Genome Biology

Alternative splicing is cell-type-specific and disease-relevant, but single-cell RNA sequencing is optimized for gene-level quantification, and how sequencing depth governs splice-junction recovery in the dominant 3'-biased chemistries has not been quantified, so experiments cannot be designed to a depth target and datasets cannot be interpreted against expected recovery. We generated deep single-cell RNA sequencing of human airway epithelial cells (12 samples, 102 to 803 million uniquely mapped reads per sample) and benchmarked junction recovery in five annotated samples against matched bulk RNA sequencing from a 190-donor cohort. Junction detection did not saturate with depth, and validity was governed by read support rather than annotation status. Recovering 90% of the observed full-depth junction complement required approximately 52,250 reads per cell, several-fold more than gene-level analysis. Quantifying splicing was more depth-limited than detecting it, and rare cell types were limited by cell number rather than depth. Against matched bulk, single-cell recall plateaued at a practical ceiling set by sample size and sparse sampling, and was independent of distance from the 3' end. These results establish a framework for designing splicing experiments in deep 3'-biased single-cell RNA sequencing, and provide an airway-epithelial dataset with matched bulk for benchmarking splicing methods.

Noncoding regulatory elements in the human genome can control gene expression through distinct mechanisms, including acting as enhancers, silencers, or topological elements. However, we lack tools to identify all such classes systematically and quantify their effects on gene expression in an endogenous genomic context. Here we develop Swap-seq, a method that combines pooled twin prime editing with flow-sorting to delete hundreds of endogenous regulatory sequences and measure their quantitative effects on the expression of a nearby gene. We apply Swap-seq to dissect the regulatory elements that control PPIF expression in THP-1 monocytes. We quantify the effects of distal and intronic enhancers, characterize CTCF insulators, and discover a silencer that regulates distal genes through histone deacetylase (HDAC)-mediated repression of neighboring enhancers. Hundreds of other genomic elements share its chromatin signatures and transcription factor binding sites, suggesting that this silencer represents a broader class of repressive elements. We leverage Swap-seq data to benchmark sequence-to-function models and show that they often capture the local activity of regulatory elements but fail to link them to distal target genes. Swap-seq thus provides a generalizable tool to discover regulatory elements across mechanistic classes, characterize their effects on gene expression, and generate large-scale datasets for evaluating sequence-to-function models.

Ribosome profiling experiments capture a snapshot of translation by sequencing ribosome-protected fragments (RPFs). Alignment to a well-curated transcriptome captures most RPF signal, but can miss ribosome activity when expressed exons are absent from the reference. Genome alignment of RPFs expands the detectable sequence space, but can expose ambiguity between homologous loci. Existing workflows often couple reference choice with the handling of ambiguously mapped reads, making it difficult to isolate the effects of either decision on downstream analyses. Here, we developed RiboFlow v2 by adding genome alignment and configurable handling of multimapping reads to the existing transcriptome workflow. Using 24 human ribosome profiling libraries and their matched RNA-seq data, we found that nucleotide-level coverage at matched coding-sequence positions and translation efficiency (TE) were broadly concordant between alignment strategies under stringent alignment filters. Even under these conditions, individual genes showed substantial differences in coverage and TE estimates. These differences included coverage in alternative exons absent from selected transcripts and ambiguous alignments between protein-coding genes and processed pseudogenes. Clustering genes by read-assignment patterns identified groups associated with pseudogene prevalence and highlighted shared exonic sequence among reference transcripts as a source of alignment discrepancies. RiboFlow v2 extends the analysis beyond selected reference transcripts and enables identification of genes whose ribosome coverage and abundance estimates depend on reference composition and how ambiguously mapped reads are handled.

Single-cell RNA sequencing (scRNA-seq) studies often profile only one or a few biological samples, limiting assessment of discovery stability under between-sample variation. We present scStable, a framework for assessing scRNA-seq analysis stability by generating synthetic samples that preserve cell-level heterogeneity from a single-sample reference while injecting bulk-derived between-sample variation. In gastric cancer data comprising 21 matched scRNA-seq and bulk RNA-seq profiles, scStable better recapitulates observed between-sample variation than ZINB-WaVE, subsampling, and Gaussian noise injection. Using synthetic sam-ples, scStable enables stability-aware analysis, including prioritization of stable differentially expressed genes and selection of clustering methods and hyperparameters.

Droplet-based and split-pool barcoding are the two dominant high-throughput single-cell RNA-seq paradigms, yet their systematic biases remain incompletely characterized. We compared seven matched 10X Genomics Chromium and Parse Biosciences Evercode WT sample pairs from five studies spanning multiple tissues and species focusing on gene counts and compositional differences rather than cell-type resolution. By separating Parse reads by priming strategy, we distinguished biases arising from priming chemistry from platform-level effects. These effects were largely independent, with platform-level differences substantially larger. Priming strategy introduced a GC-content bias within Parse, whereas the two platforms showed a consistent compositional shift across priming strategies, supported by differential capture of nuclear- and cytoplasmic-enriched transcripts. We interpret this shift as a relative enrichment of nuclear RNA in Parse. Finally, using contamination measures that do not assume a droplet-specific model, we found background RNA contamination in Parse to be comparable in magnitude to that in 10X, contrary to claims of substantially lower contamination. Together, these results indicate that expression measurements and cell-type composition are not directly interchangeable between the two platforms.

Soil ecosystems represent a vast, largely uncharacterised reservoir of microbial diversity. Metagenomic assembly has unlocked access to this resource, and recent advances in long-read sequencing have improved the recovery and quality of microbial genomes from complex samples. While genomic language models have proven highly effective at capturing biological concepts from large sequence datasets, they have predominantly been trained on reference genomes or short-read assemblies. Here, we present LOAM, a family of decoder-only genomic language models ranging from 25 to 624 million parameters and trained exclusively on Oxford Nanopore long-read environmental metagenomes. LOAM model performance scales predictably with model size and training-token budget. Despite a relatively small training sequence corpus, LOAM models outperformed comparably sized models across biological benchmarks, and achieved performance competitive with substantially larger state-of-the-art models trained on much larger datasets. Context-intervention experiments further showed that LOAM models use genomic information over several kilobases, highlighting the potential value of increased contiguity provided by long-read metagenome-assembled genomes. For probe-based benchmark tasks, we systematically evaluated representations across hidden layers and revealed that task-relevant biological information was frequently more linearly accessible from intermediate than final model layers. Finally, we observed that variation in zero-shot variant-effect prediction was strongly associated with the presence of homologous target sequences in the pre-training corpus. Together, these results establish long-read environmental metagenomes as a viable foundation for training competitive genomic language models and demonstrate the importance of both model scale and training-corpus composition for biological generalisation.

Understanding how genome sequences shape three-dimensional (3D) genome architecture is fundamental to interpreting diverse biological processes. Although previous studies have shown that sequence information can predict 3D genome architecture, these models fall short in capturing cell type-specific structures because they are trained solely on sequence inputs. The widely available Hi-C data, which contains rich structural information across biosamples, can provide complementary features to sequence data for studying cell type-specific architectures. Recently, DNA foundation models have demonstrated encouraging performance in capturing long-range genomic dependencies, holding promise for modeling chromatin interactions. However, the extremely high computational cost of running these models limits their applicability to Hi-C analysis, which requires genome-wide sequence embeddings. Here, we present Evo2HiC, a multimodal foundation model that jointly models genome sequences and structures to study cell type-specific chromatin structure. The key idea of Evo2HiC is to distill a large-scale DNA foundation model, Evo 2 (7B), into a compact encoder, while guiding the distillation with Hi-C data to preserve genomic features critical for 3D genome analysis. The model supports two types of encoders, one that operates directly on DNA sequences, and a second that additionally takes as input corresponding Hi-C data. Using the DNA-only encoder and predicting raw Hi-C contact matrices, Evo2HiC improved distance-stratified Spearman correlation by 10.9% over Orca. Moreover, by jointly embedding Hi-C and sequence information Evo2HiC achieved the best overall Pearson correlation when predicting five representative epigenomic assays. Interpretation analysis of Evo2HiC revealed its ability to identify cell type-specific sequence motifs that explain changes in epigenomic signals. Finally, we demonstrated the cross-species generalizability of Evo2HiC on 177 species from the DNA Zoo dataset for Hi-C resolution enhancement. In summary, Evo2HiC is a foundation model that integrates genome sequences and 3D chromatin structure information, substantially reduces computational cost while maintaining state-of-the-art accuracy on predicting various epigenomic signals and genome architecture, enables the identification of cell type-specific motifs, and demonstrates robust generalizability across species.

Long-read sequencing captures full-length transcripts but suffers elevated error rates, and PCR duplicates can distort quantification unless accounted for with UMIs. We present Nailpolish, a fast reference-free tool that error-corrects and deduplicates UMI-tagged long reads by generating consensus sequences via partial order alignment, while separating false duplicates arising from UMI collisions. Across eleven Oxford Nanopore Technologies and Pacific Biosciences datasets spanning bulk, single-cell, spatial, and targeted protocols, Nailpolish reduced per-read error rates and outperformed existing deduplication tools.

DNA methylation modulates transcription factor (TF) binding, and plays an important role in regulating gene expression. We developed methylTFR, an R package for quantifying TF activity from DNA methylation across TF-binding sites, providing a low-dimensional and interpretable activity profile. Across 147 human immune cell methylomes, methylTFR identified cell-type-specific regulators, including CEBP and ETS family activity in myeloid cells. In CD4+ T cells, TF activity scores resolved the naive-to-memory trajectory and highlighted AP-1 factors as important memory regulators. We applied methylTFR to sparse single-cell methylome data, and we integrated derived activity scores with gene expression and accessibility into joint factor models yielding multimodal views of lineage regulators.

Spatial transcriptomics is increasingly used to infer cellular phenotypes, molecular interactions, and functional tissue organization from RNA, although most of these biological processes are ultimately mediated by proteins. Here, we directly tested the interchangeability of spatial RNA and protein measurements using paired Xenium profiling of 21 matched targets across ~1.8 million cells from primary lung adenocarcinoma and metastatic castration-resistant prostate cancer. After accounting for marker- and dataset-specific protein background, we quantified concordance at multiple levels, including single-cell detection, expression magnitude, cell-type distribution, patient-level abundance, and spatial organization. Agreement at single-cell level was strongly marker-dependent. Canonical lineage markers generally showed concordant RNA and protein patterns, whereas several molecules, including in our study PD-1, PD-L1, and LAG-3, showed little concordance across modalities. Importantly, spatial organization remained substantially more discordant between modalities than overall abundance, indicating that spatial biomarkers are not directly transferable across modalities. Together, our results show that spatial RNA and protein measurements provide largely overlapping information for lineage markers and overall abundance, but can diverse for specific molecules and for higher-order spatial features, where the two layers provide complementary rather than equivalent readouts. These findings provide practical guidance on when spatial transcriptomic readouts can be interpreted as proxies for protein and when orthogonal protein validation is advisable.

Nucleosomes regulate DNA accessibility and thereby influence gene expression, DNA repair, and mutagenesis. Although DNA sequence strongly determines nucleosome positioning, the higher-order sequence grammar underlying nucleosome organization and its relationship to cancer-associated mutational and epigenetic processes remain poorly understood. We developed DeepND, an interpretable deep learning framework that distinguishes nucleosomal from inter-nucleosomal DNA and identifies sequence features associated with nucleosome architecture. Our interpretable framework enabled analysis of sequence patterns and identified 210 sequence motifs of 2-16 bp length that were significantly enriched or depleted in nucleosomal DNA, substantially extending the previously characterized nucleosomal motifs. These motifs revealed distinct sequence grammars within nucleosomal footprints and linker regions, including patterns associated with inward- and outward-facing minor grooves. Forty-four of these motifs showed cancer-type-specific enrichment or depletion of somatic mutations, and six motifs were associated with differential DNA methylation levels. In lung and endometrial cancers, mutations exhibited a significant ~10-bp periodic enrichment pattern across distinct minor-groove orientations, consistent with differential sensitivity of mutational and DNA repair processes to nucleosomal rotational positioning. We further identified sequence-specific methylation patterns consistent with reduced accessibility of nucleosomal DNA to DNA methyltransferases, including differential methylation at inward- and outward-facing CpG motifs. Together, these findings provide a sequence-resolved characterization of nucleosome-level organization and its link to cancer-specific mutational processes and epigenetic landscapes.

Selection of an appropriate coding exon is an important upstream step in exon deletion-based loss-of-function strategies, particularly for genes with multiple transcript isoforms. However, manual comparison of exon structures across multiple transcripts can be time-consuming and may overlook transcript-specific differences in exon boundaries and coding sequence (CDS) composition. Non3nExonFinder was developed as an interactive R/Shiny-based tool for identifying shared non-triplet coding exons across user-defined human or mouse NCBI RefSeq protein-coding transcripts. Non3nExonFinder evaluates exon/CDS annotations and coding contributions for each transcript, identifies common exons with identical genomic boundaries, and prioritizes candidates based on CDS-only status and coding-length divisibility by three. Non-triplet CDS-only exons that are internal in all selected transcripts are classified as high-confidence candidates. Independent reconstruction from the original RefSeq GTF files across 80 genes (40 human and 40 mouse), comprising 6,848 transcript-specific exon instances, showed complete concordance with Non3nExonFinder outputs, with exact recovery of all candidate sets. Non3nExonFinder provides a reproducible upstream framework for transcript-aware prioritization of exon-deletion targets.

Diploid genome assemblies are now routinely available, but most read aligners were designed for haploid references, which have long been the gold standard. When reads are aligned to a diploid assembly, the aligner sees two nearly identical alignments to either haplotype, thus reducing the mapping quality (MapQ) score to reflect this ambiguity. This can cause downstream tools to discard reads from easily mappable regions. Here, we present HipHap (HIgh- Performance HAPlotype assigner) to resolve this issue by aligning reads to each haplotype assembly separately and assigning each read to its best-supported haplotype. HipHap increased the fraction of reads mapping to the diploid genome with high MapQ, outperforming alignment to either parental haplotype. Additionally, we introduce a haplotype assignment quality score (HapQ) in HipHap to quantify confidence in the haplotype of origin of a read. HipHap is implemented in Rust, supports SAM, BAM, CRAM, and PAF formats, and is freely available at: https://github.com/jheinz27/hiphap.

Reliable classification of transposable elements (TEs) is important for genome annotation and evolutionary analysis, yet remains challenging: TE families differ widely across taxa, reference databases label them inconsistently, and many TEs are represented only by degraded copies. We introduce GLARE (Genomic Language model for Automated Repetitive Element classification), a TE classifier built by fine-tuning a genomic language model on harmonized PanTEon and Repbase reference data to predict TE superfamily across 32 superfamilies in 11 orders. Among the genomic language models compared, Nucleotide Transformer v3 provided the strongest representation. Built on this backbone, GLARE exceeded all nine classifiers in the PanTEon benchmark, reaching 0.966 weighted superfamily F1 on its 30-superfamily task. GLARE also showed stronger generalization to clades absent from training, consistently exceeding state-of-the-art classifiers on holdouts spanning three taxonomic ranks: the order Poales, the phylum Arthropoda, and the kingdom Fungi. Training on both PanTEon and Repbase reduced the performance drop seen when a single-source model was applied to the other database. Under simulated sequence decay, GLARE was the most robust classifier tested, and it performed best on independently curated degraded libraries. These findings show that genomic language models fine-tuned on harmonized multi-source references can classify TEs robustly under the taxonomic, database, and sequence-quality shifts that matter in practice. GLARE is available at https://code.jgi.doe.gov/NHaas/glare.

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

Background: Protein-coding DNA exhibits characteristic three-base periodicity. In unaligned DNA libraries, however, this signal should cancel if fragment boundaries are uniformly distributed across codon phases. Yet, we find that in multiple metagenomic libraries, base frequencies depended systematically on distance from fragment boundaries, producing a period-3 pattern. We investigated the origin of this unexpected signal. Results: The signal was observed in modern metagenomic datasets and widely across ancient-DNA datasets. We developed a new alignment-free algorithm to infer latent triplet-phase structure from nucleotide composition. The inferred shifts were strongly nonuniform at fragment boundaries, explaining the read-coordinate periodicity. We initially suspected a protocol artifact. To test an alternative, we derived a probabilistic framework linking nucleotide-dependent boundary selection to codon-phase frequencies. The framework showed that, because nucleotide composition differs among codon positions, preferential breakage at particular nucleotides biases the phases represented at fragment boundaries. Simulations confirmed this mechanism: purine-associated fragmentation produced strong period-3 coherence in a coding-dense bacterial genome, whereas no comparable genome-wide effect appeared in the human reference genome. Conclusions: The observed read-coordinate three-base periodicity does not require an intrinsic phase bias in the source DNA. Instead, sequence-dependent fragment formation or recovery can couple boundaries to codon phase. Purine-associated post-mortem fragmentation plausibly explains the widespread phenomenon in ancient-DNA libraries, although PCR- or library-specific mechanisms may contribute to it in modern datasets. Our alignment-free algorithm provides a practical method for estimating and normalizing phase shifts before searching for positional nucleotide patterns, preventing boundary-associated bias from being mistaken for an intrinsic biological pattern.

Standard spatial transcriptomics workflows rely on clustering methods to define regions of interest for downstream analyses, such as differential expression and pathway enrichment. Although researchers are presented with a single partition, most clustering methods are stochastic, meaning their output depends on random algorithmic choices, typically controlled by a seed; thus, repeating the analysis changing only the seed can yield a different result. Spatial transcriptomics data may be especially sensitive to algorithmic perturbations due to high dimensionality, sparsity, and mixed expression profiles at tissue boundaries and transitional zones; however, the extent of seed-dependent variation and its downstream impact have yet to be systematically evaluated. To bridge this gap, we present InSTability, a framework that localizes disagreement among partitions at the observation level. InSTability is built on (1) an interpretable, per-spot score that quantifies the variability of a spot's cluster assignment across partitions, and (2) an algorithm with time complexity linear in the number of spots for a specified number of runs for computing the instability score at all spots. Together, these advances enable spatially-resolved clustering instability maps for large, high-resolution datasets. Applying InSTability to repeated runs of seven established methods (BayesSpace, GraphST, Leiden, Louvain, SEDR, STAGATE, and SpiceMix) on 18 benchmark samples reveals that clustering instability varies across tissue sections, with shared spatial patterns emerging among different methods applied to the same sample. Furthermore, clustering instability correlates with poor agreement against reference annotations and propagates downstream into differential expression analyses, resulting in inconsistent gene sets across seeds. Thus, InSTability can strengthen spatial transcriptomics studies by highlighting regions where clustering is seed-sensitive and further validation is warranted. InSTability is open-source and available at https://github.com/molloy-lab/InSTability.

Background: Ancient microbial evidence informs research on past health, diet, and environments, but its scarcity and susceptibility to contamination demand care in selecting and processing data. Language models can now turn a research request into an analysis plan, but an unconstrained planner may reference files that do not exist, apply operations a format forbids, or answer a question the data cannot settle. A planner used on such evidence therefore needs to be bounded, and its decisions recoverable afterwards. Results: We developed PaleoRigor around three commitments. Agent: a language model plans the analysis. Constraint: a validator checks every file reference, format and step dependency before anything runs, execution is restricted to predefined local modules, and a planning contract directs the model to refuse out-of-scope requests with a named reason code rather than answer them. Reproducibility: each run is exported with its plan, intermediate files, checksums and decisions. In a preregistered, frozen held-out evaluation, PaleoRigor passed 23 of 24 runs (95.8%; Wilson 95% confidence interval, 79.8-99.3%): 11 of 12 supported runs and all 12 prespecified boundary decisions. An ablation that removed this planning contract from the same model passed 19 of 24 (79.2%). Repeating the frozen design with a second model configuration yielded 24 of 24 successes for PaleoRigor and 19 of 24 for the ablated arm. Pooled post hoc across both configurations, PaleoRigor passed 47 of 48 runs (97.9%), compared with 38 of 48 (79.2%) for the ablated arms. All 96 run labels were recomputed from retained records and matched. Dataset-level tests matched read counts and compressed-file sizes for six public sequencing records and documented the removal of 114 duplicate sequence entries. A retraction-associated audit illustrated how expert review limited conclusions to file-level evidence. Conclusions: For paleomicrobiome data, pairing a language-model agent with explicit constraints and a reproducible record makes the agent's proposed analysis reviewable before its outputs support a biological claim. PaleoRigor delivers checked inputs, a recorded transformation history and a named reason for each request it refuses, and runs as a local desktop application. It prepares data for specialist analysis; ancient-DNA authentication and biological interpretation remain with the paleobiologist.

Functionally linking enhancers to gene promoters at high resolution often entails analysis of 3C datasets, including viewpoint-based methods such as Capture-(Hi)C. Quantitative comparison across samples and experiments, together with visualization and integration of promoter interactions with other datasets are of critical importance. While existing statistical tools identify significant interactions and quantify underlying signal, systematically comparing the shape and extent of this signal across large sets of intervals remains challenging. Moreover, experimental artefacts, such as inefficient crosslinking, can affect downstream interpretation, but are difficult to detect systematically. Here we introduce PostChicago, an R/Bioconductor package that fills this gap. Building upon the widely used interaction-calling package Chicago, PostChicago provides a complete workflow for quality control, visualization, integration and quantitative analysis of promoter interactions between multiple samples, including annotation and quantification of genomic regions based on external datasets. PostChicago extends visual quantitative comparisons beyond individual viewpoints by computing normalized metaprofiles around large sets of intervals that visualize the distribution of interaction signal, enabling comparison between narrow and broad interaction patterns. We further demonstrate how PostChicago can identify experimental crosslinking artefacts. Together, PostChicago enables integrative, shape-sensitive analysis and quality control of promoter Capture-(Hi)C datasets at a scale that would otherwise require extensive custom scripting.

As benchmarks saturate, traditional validation methods fail to capture the probabilistic representations of genomic foundation models. This study probes self-supervised generalists (NTv3 650M, Evo2 40B, Genos 10B) alongside task-specific specialists (SpliceAI, AlphaGenome, Borzoi) using a mechanistically focused paradigm centered on expert-curated single-nucleotide substitutions within the splicing regions of the exon-rich OPA1 gene. Extending this analysis across the broader corpus of OPA1 mutations, variants of uncertain significance (VUS), and an independent 65-gene dataset of spliceogenic mutations demonstrates that generalists can match task-specific networks. Training data scale, genetic diversity, and evolutionary context are key for capturing RNA splicing syntax, whereas sparse routing may constrain performance. Mechanistically, foundation models exhibit a compute-asymmetry: the Evo2 40B model frontloads compute into massive parameter scale to resolve single-nucleotide variants within a narrow 425-bp context window, whereas the bidirectional NTv3 650M model and the Genos 10B mixture-of-experts (MoE) architecture backload compute to inference, requiring spatial logit aggregation and expansive context windows of up to 94 kb. Furthermore, classification accuracy across these foundation models benefitted from locus-specific thresholding rather than mutation-class calibration, exposing representation boundaries of the models.